{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# High-Precision Propagation\n",
"\n",
"This example demonstrates high-precision orbit propagation using satkit's advanced numerical integration capabilities. It propagates a GPS satellite orbit over a 24-hour period and validates the results against high-fidelity SP3 ephemeris data from the European Space Agency (ESA).\n",
"\n",
"## Overview\n",
"\n",
"1. **Load reference data** \u2014 read GPS satellite positions from an SP3 file of precise orbit determination results (positions only, no velocity).\n",
"2. **Rotate to inertial frame** \u2014 transform the ITRF positions into GCRF.\n",
"3. **Fit initial conditions** \u2014 use a linearized least-squares loop over the state transition matrix to recover the 6-DOF initial state, wrapped in a non-linear solve for the solar radiation pressure coefficient ($C_r A/m$) and a constant empirical acceleration in the NTW frame.\n",
"4. **Validate** \u2014 compare the propagated positions against SP3 truth and plot the residuals.\n",
"5. **Fit vs. free-run** \u2014 fit over day 1 only, then propagate through day 2 with no new measurements to see how errors grow outside the fit window.\n",
"6. **Compare integrators** \u2014 benchmark several Runge-Kutta integrators and the Gauss-Jackson 8 multistep predictor-corrector on the same problem.\n",
"\n",
"The SP3 file contains a full 24 hours of satellite positions, recorded once every 5 minutes."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Imports\n",
"import gzip\n",
"import os\n",
"import time\n",
"import urllib.request\n",
"\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import scienceplots # noqa: F401\n",
"from scipy.optimize import minimize\n",
"\n",
"import satkit as sk\n",
"\n",
"plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
"%config InlineBackend.figure_formats = ['svg']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup and SP3 File Reader\n",
"\n",
"The SP3 file format contains precise satellite positions (and optionally velocities) at regular intervals, typically generated by post-processing networks of ground stations. We read the file and extract ITRF positions for a single GPS satellite."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Function to read in the SP3 file\n",
"\n",
"\n",
"def read_sp3file(fname, satnum=20):\n",
" \"\"\"\n",
" Read SP3 file\n",
" (file containing \"true\" GPS ephemerides)\n",
" and output UTC time and position in ITRF frame\n",
" \"\"\"\n",
"\n",
" # Read in the test vectors\n",
" with open(fname, \"r\") as fd:\n",
" lines = fd.readlines()\n",
"\n",
" def line2date(lines):\n",
" for line in lines:\n",
" year = int(line[3:7])\n",
" month = int(line[8:10])\n",
" day = int(line[11:13])\n",
" hour = int(line[14:16])\n",
" minute = int(line[17:19])\n",
" sec = float(line[20:32])\n",
" yield sk.time(year, month, day, hour, minute, sec)\n",
"\n",
" def line2pos(lines):\n",
" for line in lines:\n",
" lvals = line.split()\n",
" yield np.array([float(lvals[1]), float(lvals[2]), float(lvals[3])])\n",
"\n",
" datelines = list(filter(lambda x: x[0] == \"*\", lines))\n",
" match = f\"PG{satnum:02d}\"\n",
" satlines = list(filter(lambda x: x[0:4] == match, lines))\n",
" dates = np.fromiter(line2date(datelines), sk.time)\n",
" pitrf = np.stack(np.fromiter(line2pos(satlines), list), axis=0) * 1.0e3 # type: ignore\n",
"\n",
" return (pitrf, dates)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Load SP3 Data and Seed the Initial State\n",
"\n",
"Download the day-1 SP3 file, rotate the ITRF positions into GCRF, and build a crude initial state from the first two points. The velocity seed is just a finite difference across 5 minutes \u2014 coarse, but good enough for the fit loop to converge."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Download SP3 file if not present\n",
"fname = \"./ESA0OPSFIN_20233640000_01D_05M_ORB.SP3\"\n",
"url = \"http://navigation-office.esa.int/products/gnss-products/2294/ESA0OPSFIN_20233640000_01D_05M_ORB.SP3.gz\"\n",
"\n",
"if not os.path.exists(fname):\n",
" with urllib.request.urlopen(url) as response:\n",
" with open(fname, \"wb\") as out_file:\n",
" with gzip.GzipFile(fileobj=response) as uncompressed:\n",
" out_file.write(uncompressed.read())\n",
"\n",
"# Read in the SP3 file\n",
"[pitrf, timearr] = read_sp3file(fname)\n",
"\n",
"# Rotate positions to the GCRF frame\n",
"pgcrf = np.stack(\n",
" np.fromiter(\n",
" (q * p for q, p in zip([sk.frametransform.rotation(sk.frame.ITRF, sk.frame.GCRF, t) for t in timearr], pitrf)),\n",
" list, # type: ignore\n",
" ),\n",
" axis=0,\n",
") # type: ignore\n",
"\n",
"# Crude estimation of initial velocity from the first two position states.\n",
"# The 5-minute spacing makes this rough, but the fit loop will refine it.\n",
"vgcrf = (pgcrf[1, :] - pgcrf[0, :]) / (timearr[1] - timearr[0]).seconds\n",
"\n",
"# Initial 6-DOF state (position + velocity) in GCRF\n",
"state0 = np.concatenate((pgcrf[0, :], vgcrf))\n",
"\n",
"# Propagator settings shared across fits.\n",
"# Available integrators: sk.integrator.rkv98 (default), rkv98_nointerp, rkv87, rkv65, rkts54\n",
"# Available gravity models: sk.gravmodel.egm96 (default), jgm3, jgm2, itugrace16\n",
"settings = sk.propsettings(\n",
" abs_error=1e-12,\n",
" rel_error=1e-12,\n",
" gravity_degree=10,\n",
")\n",
"# Only compute sun, moon positions and earth rotation vectors once for all propagations\n",
"settings.precompute_terms(timearr[0], timearr[-1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit Initial State and Validate\n",
"\n",
"We fit the orbit in two nested loops:\n",
"\n",
"- **Inner loop (linearized least squares)** \u2014 iteratively solves for the 6-element initial state correction $\\Delta s_0$ using the position block of the state transition matrix $\\Phi_p(t; s_0)$:\n",
"\n",
"$$\\min_{\\Delta s_0} \\; \\sum_t \\| p_\\text{gps}(t) - \\hat{p}(t; s_0) - \\Phi_p(t; s_0)\\,\\Delta s_0 \\|^2$$\n",
"\n",
" where $p_\\text{gps}(t)$ is the SP3 position in GCRF and $\\hat{p}(t; s_0)$ is the propagated position from initial state $s_0$. Converges in about 5 iterations.\n",
"\n",
"- **Outer loop (Nelder-Mead)** \u2014 non-linearly searches for parameters the STM doesn't cover: the solar radiation pressure coefficient $C_r A/m$ and a **constant empirical acceleration** in the velocity-aligned NTW frame $[a_N, a_T, a_W]$.\n",
"\n",
"The empirical acceleration is the interesting physical piece. It is **not** a real thruster \u2014 it is a catch-all for any constant bias left over after the gravity / third-body / SRP model is applied. Such biases come from things like mismodelled SRP box-wing geometry, antenna thrust, thermal re-radiation, truncation of the spherical-harmonic gravity expansion, or errors in Earth-orientation data. Because we use the NTW frame, the tangential component $a_T$ acts purely along the velocity vector and therefore maps one-to-one to orbital energy (and hence semi-major-axis drift) \u2014 it is the most effective single parameter for absorbing along-track residual growth, which is the dominant error mode in GNSS orbit determination. In operational orbit determination these are called \"empirical accelerations\" and fitting them is standard practice (e.g. the CODE 5-parameter SRP model, JPL's GPS orbits, etc.)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Fit helpers: satproperties builder and linearized LSQ loop\n",
"\n",
"\n",
"def make_satprops(craoverm, ntw_accel):\n",
" \"\"\"Build a satproperties with a constant NTW-frame thrust acceleration.\n",
"\n",
" NTW axes (Vallado \u00a73.3): T = velocity unit vector, W = orbit-normal\n",
" (r \u00d7 v), N = W \u00d7 T (in-plane, normal to velocity). Strict along-velocity\n",
" semantics even for eccentric orbits \u2014 the T component acts purely on\n",
" orbital energy (semi-major axis).\n",
" \"\"\"\n",
" thrusts = [\n",
" sk.thrust.constant(\n",
" list(ntw_accel), timearr[0], timearr[-1], frame=sk.frame.NTW\n",
" )\n",
" ]\n",
" return sk.satproperties(craoverm=craoverm, thrusts=thrusts)\n",
"\n",
"\n",
"def linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp, iters=5):\n",
" \"\"\"\n",
" Linearized least squares fit of initial state to SP3 truth data, holding\n",
" CrAoverM and any thrust terms fixed.\n",
" \"\"\"\n",
" state0_s = state0.copy()\n",
" for idx in range(iters):\n",
" # Propagate state and state transition matrix over times of interest\n",
" res = sk.propagate(\n",
" state0_s,\n",
" timearr[0],\n",
" timearr[-1],\n",
" output_phi=True,\n",
" propsettings=settings,\n",
" satproperties=sp,\n",
" )\n",
"\n",
" # Get state and state transition matrix at times of GPS truth data\n",
" statearr, phiarr = zip(*[res.interp(t, output_phi=True) for t in timearr])\n",
" phiarr = np.array(phiarr)\n",
" statearr = np.array(statearr)\n",
"\n",
" # Linearized least squares solve for state0 update\n",
" H = np.sum([p[0:3, :].T @ p[0:3, :] for p in phiarr], axis=0)\n",
" b = np.sum(\n",
" [\n",
" p[0:3, :].T @ (pgcrf[i, :] - statearr[i, 0:3]).T\n",
" for i, p in enumerate(phiarr)\n",
" ],\n",
" axis=0,\n",
" )\n",
" dstate0 = np.linalg.solve(H, b)\n",
" state0_s = state0_s + dstate0\n",
"\n",
" perr = np.zeros((len(timearr), 3))\n",
" for i in range(len(timearr)):\n",
" perr[i, :] = res.interp(timearr[i])[0:3] - pgcrf[i, :]\n",
"\n",
" return state0_s, res, perr"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Outer loop: fit $C_r A/m$ and empirical NTW acceleration\n",
"\n",
"Nelder-Mead searches over four scalars: $C_r A/m$ and $[a_N, a_T, a_W]$. At each outer step the inner linearized LSQ re-fits the 6-element initial state (the STM only covers position/velocity, so these acceleration parameters must be handled by the outer non-linear loop). The NTW accel components are scaled by `ACCEL_SCALE` so all parameters are $O(0.01)$ and Nelder-Mead behaves."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Nelder-Mead outer optimization over CrA/m and constant NTW empirical acceleration\n",
"ACCEL_SCALE = 1e-9 # m/s^2 per unit of scaled parameter\n",
"\n",
"\n",
"def minfunc_full(params, state0, timearr, pgcrf, settings):\n",
" craoverm = params[0]\n",
" ntw_accel = np.array(params[1:4]) * ACCEL_SCALE\n",
" if craoverm < 0 or craoverm > 1:\n",
" return 1e30\n",
" sp = make_satprops(craoverm, ntw_accel)\n",
" _, _, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp)\n",
" return float(np.sum(perr**2))\n",
"\n",
"\n",
"x0 = np.array([0.01, 0.0, 0.0, 0.0])\n",
"r = minimize(\n",
" minfunc_full,\n",
" x0,\n",
" args=(state0, timearr, pgcrf, settings),\n",
" method=\"Nelder-Mead\",\n",
" options={\"xatol\": 1e-5, \"fatol\": 1e-4, \"maxiter\": 400},\n",
")\n",
"\n",
"craoverm_fit = r.x[0]\n",
"ntw_accel_fit = np.array(r.x[1:4]) * ACCEL_SCALE\n",
"\n",
"# Final least-squares fit with optimized CrA/m and NTW empirical acceleration\n",
"sp = make_satprops(craoverm_fit, ntw_accel_fit)\n",
"state0, res, perr = linearized_least_squares_fit(state0, timearr, pgcrf, settings, sp)\n",
"\n",
"print(f\"Satellite radiation pressure susceptibility, Cr A over M: {craoverm_fit:.4f} m^2/kg\")\n",
"print(\n",
" f\"Fitted empirical NTW acceleration (m/s^2): \"\n",
" f\"N={ntw_accel_fit[0]: .3e}, T={ntw_accel_fit[1]: .3e}, W={ntw_accel_fit[2]: .3e}\"\n",
")\n",
"print(\n",
" f\"Mean position error of fit over 24 hours: {np.mean(np.linalg.norm(perr, axis=1)):.3f} meters\"\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Validate: residuals vs SP3\n",
"\n",
"Per-component position errors relative to the SP3 truth over the full 24-hour fit window."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Plot position error\n",
"fig, ax = plt.subplots()\n",
"dates = [t.datetime() for t in timearr]\n",
"ax.plot(dates, perr[:, 0], label=\"$X$\")\n",
"ax.plot(dates, perr[:, 1], label=\"$Y$\")\n",
"ax.plot(dates, perr[:, 2], label=\"$Z$\")\n",
"ax.set_xlabel(\"Time\")\n",
"ax.set_ylabel(\"Position Error (m)\")\n",
"ax.set_title(\"Propagation Error vs SP3 Truth for GPS Satellite\")\n",
"ax.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.25), ncol=3)\n",
"fig.autofmt_xdate()\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fit vs. Free-Run: How Errors Grow\n",
"\n",
"The least-squares fit above stayed at the meter level over its 24-hour window. But what happens if we hold the fitted initial state fixed and propagate *beyond* the fit interval?\n",
"\n",
"Below we download a second day of SP3 truth data and propagate the already-fitted `state0` across both days. The day-1 window is the fit (errors at the meter level); day 2 is pure prediction with no new measurements, and shows how quickly errors grow once we leave the fit interval."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Download a second SP3 file (day after the first) and concatenate to get 2 days of truth data.\n",
"# Day 1 is 2023-364 (GPS week 2294). Day 2 is 2023-365 = Dec 31 2023 (Sunday), which starts GPS week 2295.\n",
"fname2 = \"./ESA0OPSFIN_20233650000_01D_05M_ORB.SP3\"\n",
"url2 = \"http://navigation-office.esa.int/products/gnss-products/2295/ESA0OPSFIN_20233650000_01D_05M_ORB.SP3.gz\"\n",
"\n",
"if not os.path.exists(fname2):\n",
" with urllib.request.urlopen(url2) as response:\n",
" with open(fname2, \"wb\") as out_file:\n",
" with gzip.GzipFile(fileobj=response) as uncompressed:\n",
" out_file.write(uncompressed.read())\n",
"\n",
"# Read day 2 and rotate to GCRF\n",
"pitrf2, timearr2 = read_sp3file(fname2)\n",
"pgcrf2 = np.stack(\n",
" np.fromiter(\n",
" (q * p for q, p in zip([sk.frametransform.rotation(sk.frame.ITRF, sk.frame.GCRF, t) for t in timearr2], pitrf2)),\n",
" list, # type: ignore\n",
" ),\n",
" axis=0,\n",
") # type: ignore\n",
"\n",
"# Concatenate day 1 + day 2 for full 2-day truth. The last point of day 1 and\n",
"# first point of day 2 are typically the same epoch, so drop the duplicate.\n",
"if (timearr2[0] - timearr[-1]).seconds == 0:\n",
" timearr_full = np.concatenate((timearr, timearr2[1:]))\n",
" pgcrf_full = np.concatenate((pgcrf, pgcrf2[1:, :]), axis=0)\n",
"else:\n",
" timearr_full = np.concatenate((timearr, timearr2))\n",
" pgcrf_full = np.concatenate((pgcrf, pgcrf2), axis=0)\n",
"\n",
"# Propagation settings covering both days\n",
"settings2 = sk.propsettings(\n",
" abs_error=1e-12,\n",
" rel_error=1e-12,\n",
" gravity_degree=10,\n",
")\n",
"settings2.precompute_terms(timearr_full[0], timearr_full[-1])\n",
"\n",
"# Propagate the already-fitted state0 (and sp) across the full 2-day span.\n",
"# Day 1 is the fit window; day 2 is pure prediction \u2014 no new measurements used.\n",
"res_full = sk.propagate(\n",
" state0,\n",
" timearr_full[0],\n",
" timearr_full[-1],\n",
" propsettings=settings2,\n",
" satproperties=sp,\n",
")\n",
"\n",
"perr_full = np.zeros((len(timearr_full), 3))\n",
"for i, t in enumerate(timearr_full):\n",
" perr_full[i, :] = res_full.interp(t)[0:3] - pgcrf_full[i, :]\n",
"\n",
"perr_full_norm = np.linalg.norm(perr_full, axis=1)\n",
"\n",
"# Split day 1 (fit window) vs day 2 (free run) for reporting\n",
"fit_end = timearr[-1]\n",
"is_fit = np.array([(t - fit_end).seconds <= 0 for t in timearr_full])\n",
"print(f\"Mean |err| over day 1 (fit window): {np.mean(perr_full_norm[is_fit]):8.3f} m\")\n",
"print(f\"Max |err| over day 1 (fit window): {np.max(perr_full_norm[is_fit]):8.3f} m\")\n",
"print(f\"Mean |err| over day 2 (free run): {np.mean(perr_full_norm[~is_fit]):8.3f} m\")\n",
"print(f\"Max |err| over day 2 (free run): {np.max(perr_full_norm[~is_fit]):8.3f} m\")\n",
"\n",
"# Plot: components + error magnitude, shading the fit window vs free-run region\n",
"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 5), sharex=True)\n",
"dates_full = [t.datetime() for t in timearr_full]\n",
"\n",
"ax1.plot(dates_full, perr_full[:, 0], label=\"$X$\")\n",
"ax1.plot(dates_full, perr_full[:, 1], label=\"$Y$\")\n",
"ax1.plot(dates_full, perr_full[:, 2], label=\"$Z$\")\n",
"ax1.axvline(fit_end.datetime(), color=\"k\", linestyle=\"--\", linewidth=1)\n",
"ax1.axvspan(\n",
" dates_full[0], fit_end.datetime(), alpha=0.1, color=\"tab:green\", label=\"Fit window\"\n",
")\n",
"ax1.axvspan(\n",
" fit_end.datetime(), dates_full[-1], alpha=0.1, color=\"tab:red\", label=\"Free run\"\n",
")\n",
"ax1.set_ylabel(\"Position Error (m)\")\n",
"ax1.set_title(\"Fitted-State Propagation Error vs SP3 Truth (Day 1 fit, Day 2 free-run)\")\n",
"ax1.legend(loc=\"upper center\", bbox_to_anchor=(0.5, -0.12), ncol=5, fontsize=8)\n",
"\n",
"ax2.semilogy(dates_full, perr_full_norm, color=\"tab:purple\")\n",
"ax2.axvline(fit_end.datetime(), color=\"k\", linestyle=\"--\", linewidth=1)\n",
"ax2.axvspan(dates_full[0], fit_end.datetime(), alpha=0.1, color=\"tab:green\")\n",
"ax2.axvspan(fit_end.datetime(), dates_full[-1], alpha=0.1, color=\"tab:red\")\n",
"ax2.set_xlabel(\"Time\")\n",
"ax2.set_ylabel(\"|Error| (m, log)\")\n",
"\n",
"fig.autofmt_xdate()\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Comparing Integrators\n",
"\n",
"The propagator supports several Runge-Kutta integrators of different orders, plus the Gauss-Jackson 8 multistep predictor-corrector. Higher-order Runge-Kutta methods can take larger steps for the same accuracy, so they often require fewer total function evaluations despite more stages per step. Gauss-Jackson 8 is specialised for smooth 2nd-order orbit ODEs and typically needs 3\u201310\u00d7 fewer force evaluations than the adaptive methods at comparable accuracy \u2014 at the cost of a user-chosen fixed step size and no state-transition-matrix support.\n",
"\n",
"Below we compare performance and accuracy of these integrators on the same problem, using the fitted initial state from above."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Integrators to compare (from lowest to highest order). Gauss-Jackson 8 is a\n",
"# fixed-step multistep predictor-corrector \u2014 it uses gj_step_seconds from\n",
"# propsettings rather than abs_error / rel_error. 120 s is a reasonable step\n",
"# for GPS (MEO) altitudes.\n",
"integrators = [\n",
" (\"rkts54 - Tsitouras 5(4)\", sk.integrator.rkts54),\n",
" (\"rkv65 - Verner 6(5)\", sk.integrator.rkv65),\n",
" (\"rkv87 - Verner 8(7)\", sk.integrator.rkv87),\n",
" (\"rkv98 - Verner 9(8)\", sk.integrator.rkv98),\n",
" (\"gauss_jackson8 - GJ8 (120 s)\", sk.integrator.gauss_jackson8),\n",
"]\n",
"\n",
"print(\n",
" f\"{'Integrator':<35} {'Steps':>6} {'Evals':>6} {'Time (ms)':>10} {'Max Err (m)':>12}\"\n",
")\n",
"print(\"-\" * 75)\n",
"\n",
"for name, integ in integrators:\n",
" s = sk.propsettings(\n",
" abs_error=1e-12,\n",
" rel_error=1e-12,\n",
" gravity_degree=10,\n",
" integrator=integ,\n",
" gj_step_seconds=120.0, # only used by gauss_jackson8\n",
" )\n",
" s.precompute_terms(timearr[0], timearr[-1])\n",
"\n",
" t0 = time.perf_counter()\n",
" result = sk.propagate(\n",
" state0,\n",
" timearr[0],\n",
" timearr[-1],\n",
" propsettings=s,\n",
" satproperties=sp,\n",
" )\n",
" elapsed = (time.perf_counter() - t0) * 1e3\n",
"\n",
" # Compute max position error vs SP3 truth\n",
" max_err = max(\n",
" np.linalg.norm(result.interp(t)[0:3] - pgcrf[i, :])\n",
" for i, t in enumerate(timearr)\n",
" )\n",
"\n",
" stats = result.stats\n",
" print(\n",
" f\"{name:<35} {stats.num_accept:>6} {stats.num_eval:>6} {elapsed:>10.1f} {max_err:>12.3f}\"\n",
" )"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}